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1.  Introduction 


Weak  intermolecular  interactions  play  a  particularly  important  role  in  the  computer  simulation 
of  liquids  [1].  An  excellent  example  of  a  process  which  displays  this  characteristic  is  supercritical 
fluid  (SCF)  C02  extraction  [2].  When  SCFs  are  compressed  to  liquid-like  densities,  their  solvent 
strength  dramatically  increases.  A  closed  system  can  then  be  built  to  extract  materials  of  interest 
from  a  more  general  mixture.  Carbon  dioxide  (C02)  proves  to  be  an  excellent  choice  for  this  process 
because  of  both  its  nondestructive  character  and  because  it  is  environmentally  benign.  Some 
industrial  processes,  such  as  caffeine  extraction,  already  make  profitable  use  of  this  procedure. 

The  Department  of  Defense  (DOD)  community  sets  a  priority  on  developing  an  environmentally 
beneficial  and  cost-effective  way  to  recycle  solid  energetic  materials  which  have  reached  the  end  of 
their  rated  lifetime.  Significant  environmental  and  economic  advantages  could  result  from  an 
industrial  scale,  closed-system  recycling  procedure  based  on  supercritical  C02  for  this  particular 
application.  Unfortunately,  certain  components  in  composite  propellants  are  not  sufficiently  soluble 
in  pure  SCF  C02  to  make  this  extraction  process  viable.  The  solubility  characteristics  of  these 
components  can  be  enhanced  with  the  addition  of  so-called  modifier  molecules.  These  typically 
polar  molecules  increase  the  solubility  strength  of  the  SCF,  but  little  is  known  about  the  detailed 
molecular  interactions  accounting  for  the  increased  solubility.  The  first  step  towards  simulating  the 
entire  system  is  knowledge  of  accurate  intermolecular  potentials  for  all  dimer  interactions  in  the 
system — the  solvent-solute,  solvent-modifier,  modifier-solute,  and  each  with  itself  (e.g., 
solvent-solvent).  Methyl  cyanide  (CH3CN)  has  been  shown  to  be  an  effective  polar  modifier  for 
enhancing  the  dissolution  of  one  of  the  important  solid  energetic  materials, 
cyclotrimethylenetrinitramine  (RDX),  in  SCF  C02  [3] .  This  work  will  focus  on  mapping  the  detailed 
potential  energy  surface  (PES)  of  C02  interacting  with  CH3CN. 

Symmetry- Adapted  Perturbation  Theory  (SAPT)  [4, 5]  is  a  natural  choice  to  find  the  interaction 
energy  of  the  CH3CN-C02  system  or  any  two  closed-shell  weakly  interacting  atomic  or  molecular 
systems.  SAPT  directly  and  naturally  separates  the  interaction  energy  into  four  physically 
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interpretable  components:  (1)  electrostatics,  (2)  exchange,  (3)  dispersion,  and  (4)  induction.  Each 
component  has  distinct  radial  and  angular  dependence  for  each  system  and  can  be  fitted  to  an 
analytical  form  independently  of  the  other  components.  This  can  lead  to  significant  physical  insight 
about  the  interaction  in  contrast  to  the  currently  more  popular  supermolecular  (SM)  method  which 
returns  only  a  single  number.  SAPT  has  been  used  to  successfully  investigate  a  variety  of  systems 
including  Ar-H2  [6], He-HF [7], He-CO  [8], Ar-HF [9], He-C2H2  [10],  H2-CO  [11],  and (H20)2  [12]. 

Section  2  introduces  definitions  necessary  for  analyzing  SAPT  results.  Section  3  describes  the 
computational  details.  Section  4  investigates  some  representative  cuts  of  the  PES  for  the 
CH3CN-C02  system  and  specifies  the  coarse  grid  used  for  the  majority  of  the  geometrical 
configurations  investigated.  Section  5  presents  conclusions. 

2.  Method 

Jeziorski,  Moszynski,  and  Szalewicz  [4  ]  and  Szalewicz  and  Jeziorski  [5]  present  recent  reviews 
of  SAPT  and  provide  an  excellent  overview  of  the  method.  Further  details  on  the  explicit  derivation 
of  the  theory  and  implementation  can  be  found  in  section  6  [13-20].  We  present  only  a  necessary 
amount  of  notation  to  interpret  the  results  of  the  present  work. 

The  dimer  Hamiltonian  is  decomposed  by  SAPT  into  three  general  parts.  The  first  two,  the  Fock 
operator  F,  and  the  Moller-Plesset-type  intramonomer  correlation  operator  W,  have  separate 
contributions  from  both  systems  A  and  B  and  are  written  as  F  =  FA  +  FB  and  W  =  WA  +  WB , 
respectively.  The  third  part  of  the  Hamiltonian  is  the  intermolecular  interaction  operator  V  which 
mediates  interactions  between  the  two  systems.  The  total  Hamiltonian  is  then  written  as 
H  =F  +  V  +  W.  The  wave  function  used  with  .this  Hamiltonian  is  the  product  of  the  system  A  and 
B  wave  functions.  This  product  wave  function  does  not  obey  the  Pauli  principle.  The  correct 
permutational  symmetry  of  the  electrons  between  systems  is  imposed  on  the  product  wave  functions 
using  an  antisymmetrizer  described  in  more  detail  in  Jeziorski,  Moszynski,  and  Szalewicz  [4  ]  and 
Szalewicz  and  Jeziorski  [5]. 
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The  intermolecular  interaction  energy  £jnt  within  the  SAPT  framework  can  then  be  expanded 
in  powers  of  the  intermolecular  interaction  operator  V  as 
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where  the  first  two  terms  on  the  right-hand  side  of  equation  (1)  can  be  interpreted  as  the  classical 
electrostatic  (coulomb)  and  exchange  energies,  respectively.  The  exchange  components  are  the 
result  of  the  antisymmetrization  previously  mentioned.  They  can  also  be  viewed  as  an  effect  of 
resonance  tunneling  of  electrons  between  the  interacting  systems. 

The  second-order  terms  in  equation  (1)  naturally  separate  into  dispersion  and  induction 
components  as 
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and  analogously  for  the  exchange  component 

p(2)  =  Fm  .  p(2) 

^  exch  ^  exch-ind  ^  exch -disp* 


(3) 


The  dispersion  energy  is  a  result  of  the  interactions  of  the  two  monomers’  instantaneous  electric 
moments.  The  induction  energy  describes  the  interactions  of  the  permanent  and  induced  multipole 
moments  of  the  two  monomers.  The  second-order  exchange-dispersion  and  exchange-induction 
energies  result  from  electron  tunneling  between  systems  related  to  the  dispersion  and  induction 
components  of  the  wave  function. 

Equation  ( 1 )  implicitly  indicates  the  inclusion  of  full  intramonomer  electron  correlation.  Within 
the  SAPT  framework,  this  is  only  currently  possible  in  the  case  of  four-electron  systems  [21, 22]. 
To  describe  the  intramonomer  electron  correlation  for  larger  systems,  we  also  perturbationally 
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expand  each  of  the  components  in  equation  (1)  in  powers  of  W.  For  example,  the  first-order 
polarization  energy  is  now  expanded  in  a  double  perturbation  series  as 


17(1)  _  XT'  £»(!&) 
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where  k  indicates  the  order  in  W.  It  is  convenient  to  split  expansions  like  equation  (4)  into  terms 
which  include  and  neglect  intramonomer  correlation.  This  is  written  explicitly  for  the  first-order 
polarization  energy  as 
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where  the  second  term  sums  all  terms  of  order  one  and  above  in  W  in  equation  (4).  The  sum  of  these 
terms  through  the  Ath  order  in  W  will  be  indicated  by  the  notation  e^t  (A).  Similar  definitions  are 

assumed  for  the  other  components  as  well. 

The  first-order  polarization  and  second-order  induction  components  are  calculated  with  the 
inclusion  of  the  coupled  Hartree-Fock  (HF)-type  response  of  a  perturbed  system.  Components 

computed  in  this  manner  will  be  indicated  with  the  subscript  “resp”  as  in  E  (i^resp  •  The 
intramonomer  correlation  effects  in  induction  interactions  will  be  approximated  by  tE(^,  the 
so-called  “true”  correlation  contribution  which  collects  those  parts  of  the  E^  correction  that  are 
not  included  in  E ? Note  that  =  E^m p. 


The  correlated  S  APT  portion  of  the  interaction  energy  then  includes 
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The  ^exch-md  component  which  partially  quenches  the  corresponding  induction  component  is  not 

currently  coded.  We  have  estimated  it  by  scaling  E^h-ind  by  the  ratio  of  the  correlated  to 
uncorrelated  induction  components  by 
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There  exists  the  following  relation  between  the  SM  HF  interaction  energy  and  the  SAPT 
expansion  [23, 24]: 


=  E-dO)  +  E-00)  .  p  (20)  r  (20) 

int  "  elst  "  exch  ^  ind, resp  "  exch-ind, resp 


6®, 
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where  6^  indicates  the  sum  of  higher  order  induction  and  exchange  terms.  The  first-  and 
second-order  terms  in  equation  (8)  are  those  calculated  in  the  current  implementation  of  SAPT.  6^ 
is  defined  as  the  difference  between  the  sum  of  these  SAPT  terms  and  the  supermolecule  HF  energy,  E  ^ . 

In  order  to  include  some  of  the  higher  order  induction  terms  currently  not  available  in  SAPT,  we  use 
a  hybrid  method  which  includes  the  SM  HF  energy  and  the  correlated  portion  of  the  SAPT 
components  indicated  in  equation  (6).  The  total  interaction  energy  in  the  present  work  will  then  be 
approximated  by  combining  equation  (6)  and  equation  (8)  to  give 


=  E 2 


+  E 


corr 
int  • 


(9) 


For  more  discussion  about  this  relationship  and  further  references,  see  Szalewicz  and  Jeziorski  [5] . 

tip 

The  Boys-Bemardi  [25]  counterpoise  (CP)  scheme  is  always  used  to  compute  E  -mt  and  other  SM 
quantities  of  interest  in  order  to  eliminate  basis  set  superposition  error  (BSSE)  [5, 26]. 
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3.  Computational  Details 


We  used  Dunning’s  correlation-consistent  basis  augmented  with  diffuse  functions  labeled 
aug-cc-PVDZ  [27-29]  as  a  starting  point  for  all  calculations  with  modifications.  The  CH3CN 
monomer  geometry  was  determined  by  a  QCISD  [30-32]  calculation  with  the  full  inclusion  of  inner 
shell  electrons  optimizing  the  monomer’s  total  energy.  The  nuclear  coordinates  for  this  monomer 
are  presented  in  Table  1.  For  carbon  dioxide,  a  carbon-oxygen  distance  of  1.162047  A  was  taken 
from  Sadleg,  Szczesniak,  and  Chalasinski  [33].  Both  monomer  geometries  were  then  fixed  for  all 
further  study.  Gaussian  94  [34]  and  Atmol  [35]  were  both  used  to  perform  the  necessary  SCF 
calculations.  Both  programs  are  interfaced  to  the  SAPT  suite  of  codes  [36]. 


Table  1.  Nuclear  Coordinates  in  A  for  the  CH3CN  Monomer  Geometry.  Each  HCCN  Angle 
Is  109.731504°.  The  Center  of  Mass  for  the  System  is  0.68927  A  From  the  Inner 
Carbon  Between  the  Two  Carbon  Atoms. 


Atom 

X 

Z 

C 

0.0 

0.0 

0.0 

C 

1.477947 

0.0 

0.0 

N 

-1.170898 

0.0 

0.0 

H 

1.849295 

0.0 

1.0353445 

H 

1.849295 

-0.896636 

-0.517673 

H 

1.849295 

-0.896636 

-0.517673 

Optimizations  of  the  full  dimer  energy  in  the  full  dimer  basis  set  with  fixed  internal  monomer 
geometries  were  then  performed  at  the  MP2  level  of  theory  again  with  full  inclusion  of  inner  core 
electrons.  Four  local  minimum  geometries  were  located  with  this  procedure  and  will  be  designated 
Gn,  n  =  1, 2, 3,  or  4.  These  geometries,  shown  in  Figure  1,  serve  as  starting  points  for  investigating 
interesting  portions  of  the  PES.  Since  a  CP-corrected  optimization  procedure  was  not  used  for  the 
minimizations,  BSSE  may  affect  the  positions  of  the  local  minima.  No  attempt  to  quantify  this  error 
was  made  since  the  minimum  geometries  were  used  only  for  the  purpose  previously  mentioned. 
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These  Four  Minima  at  -2.90  kcal/mol. 


Subsequent  discussion  of  these  geometries  refers  to  the  SAPT  computations  which  are  not  biased 
by  BSSE. 

The  coordinates  of  the  MP2  local  minimum  geometries  are  provided  in  Table  2.  A  description 
of  the  coordinate  system  is  given  in  the  next  paragraph.  In  configuration  Gl,  the  C02  axis  is  in  a 
slipped,  nearly  parallel  position  (away  from  the  CH3  group)  with  respect  to  the  C3  axis  of  the  CH3CN 
molecule  (see  Figure  1).  In  configuration  G2,  the  C02  is  oriented  along  and  nearly  perpendicular 
to  the  major  axis  of  CH3CN  toward  the  CH3  group.  Configuration  G4  has  the  C02  in  a  similar 
location,  but  the  C02  major  axis  is  nearly  aligned  with  the  CH3CN  axis.  Finally,  configuration  G3 
is  oriented  in  the  same  way  as  G2  except  toward  the  nitrogen  atom. 


Table  2.  The  Coordinates  for  the  Local  Minimum  Geometries  Gl,  G2,  G3,  and  G4  (The  Units 
for  Distance  and  Angles  Are  A  and  Degrees,  Respectively.) 


Geometry 

R 

P. 

Yi 

P2 

«2 

Gl 

3.327156 

107.746148 

119.852998 

116.669519 

359.978866 

G2 

4.652572 

7.954955 

59.997536 

95.612989 

359.870681 

G3 

4.276012 

179.643584 

60.183120 

89.771324 

170.272920 

G4 

5.572889 

0.012626 

61.490949 

179.921826 

1.518053 

The  dimer  configuration  has  been  specified  by  coordinates  consisting  of  a  separation  distance 
R  and  two  sets  of  Euler  angles  as  given  in  Brink  and  Satchler  [37].  A  pictorial  representation  of 
these  coordinates  as  applied  to  this  system  can  be  seen  in  Figure  2.  R  is  defined  as  the  length  of  the 
vector  connecting  the  centers  of  mass  between  the  two  monomers.  Each  center  of  mass  is  located 
at  the  origin  of  a  system  of  Cartesian  coordinate  axes,  and  these  two  sets  of  axes  remain  fixed  and 
parallel  to  the  space-fixed  coordinate  system.  The  vector  R  coincides  with  the  Z  axis  and  points  in 
the  positive  Z  direction.  One  set  of  Euler  angles  is  assigned  to  each  of  the  molecules  and  is  defined 
with  respect  to  the  Cartesian  coordinate  axes  associated  with  that  molecule.  The  angles  are  given 
by  the  variables  (a2.  Pi,  Yi)  for  the  orientation  of  CH3CN,  and  (a2,  p2,  y2)  for  C02. 
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The  dimer  configurations  are  obtained  from  the  independent  Euler  angles  by  rotating  CH3CN 
about  y,,  followed  by  rotating  each  monomer  through  its  respective  pt  angle,  and  finally  rotating  the 
C02  molecular  axis  through  its  02  angle.  yY  is  the  angle  of  rotation  of  CH3CN  about  its  C3  axis 
(hereafter  referred  to  as  its  “molecular  axis”).  This  can  be  further  defined  as  a  rotation  about  the 
CCN  molecular  axis  between  a  half  plane  and  the  stationary  XZ  plane.  The  half-plane  is  formed  by 
the  CH3CN  molecular  axis  and  one  of  the  CH  bonds,  y!  is  taken  to  be  zero  when  the  half-plane 
containing  the  CH  bond  coincides  with  the  XZ  plane  and  lies  in  the  positive  X  hemisphere  (see 
Figure  2).  yl  is  then  allowed  to  vary  between  0  and  120°,  with  the  sense  of  rotation  always  being 
done  such  that  the  CH  bond  lies  in  the  positive  Y  hemisphere.  With  this  definition  of  yl5  one  can 
take  advantage  of  the  C3l)  symmetry  of  CH3CN.  y2  would  be  the  angle  of  rotation  for  C02  about  its 
molecular  axis,  but  due  to  its  cylindrical  symmetry,  it  is  undefined  and  allows  us  to  reduce  the 
coordinates  from  six  independent  coordinates  to  five. 

|3  j  is  the  simple  angle  between  the  positive  Z  axis  and  the  vector  drawn  from  the  COM  of  CH3CN 
to  the  methyl  carbon.  Likewise,  p2  is  the  angle  between  one  of  the  CO  bonds  and  the  part  of  the  Z 
axis  where  Z  >R.  Both  P!  and  P2  are  allowed  to  vary  from  0  to  180°.  02  is  an  angle  of  rotation 
between  two  half  planes,  both  of  which  have  their  edge  coincidental  with  the  Z  axis.  The  first  half 
plane,  which  is  taken  to  remain  stationary,  is  formed  by  the  Z  axis  and  the  positive  X  axis.  The 
second  half  plane  (i.e.,  the  rotating  plane)  is  formed  by  the  Z  axis  and  the  CO  bond  involved  in  the 
definition  of  p2.  02  can  vary  from  0  to  360°,  but  is  zero  when  the  half  plane  containing  the  CO  is 
coincidental  with  the  XZ  plane  and  pointing  in  the  +X  direction.  The  sense  of  rotation  is  clockwise 
when  viewed  down  the  Z  axis  from  -Z  to  +Z.  The  final  Euler  angle,  al5  would  rotate  the  molecular 
axis  for  CH3CN  out  of  the  XZ  plane.  This  corresponds  to  the  rotation  of  the  entire  dimer  system, 
in  a  fixed  configuration,  about  the  Z  axis.  Thus,  ax  is  not  involved  in  the  definition  of  the  relative 
positions  of  the  two  molecules,  and  thus  can  be  fixed  to  zero.  The  results  of  this  are  that  the  CCN 
molecular  axis  always  lies  in  the  XZ  plane,  and  the  Pt  angle  always  lies  in  the  positive  X  direction. 

There  are  a  few  symmetries  within  this  system  which  can  be  exploited  to  reduce  the  total  number 
of  grid  points  that  need  to  be  considered.  The  linearity  of  the  C02  molecule  eliminates  the  need  to 
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Coincidental  With  the  Z  Axis  and  Connects  the  Centers  of  Mass  of  the  Two  Dimers.  Associated  With  Each  Monomer 
Is  a  Set  of  Euler  Angles  as  Given  in  Brink  and  Satchler  [37].  The  Orientation  of  CH3CN  Is  Given  by  the  Angles  a„  P„ 
and  Yi»  Where  ax  Is  Held  Fixed  at  0,  and  the  Orientation  of  C02  Is  Given  by  a2  and  P2.  Due  to  the  Cylindrical  Symmetry 
of  C02,  y2  Is  Undefined,  Which  Reduces  the  Number  of  Independent  Variables  to  5.  (See  the  Text  for  Further 


consider  the  sixth  coordinate  necessary  for  describing  general  rigid  two-body  interactions. 
Describing  £'int  as  a  potential  function  of  five  coordinates  by 


£*  -  (10) 

the  symmetries  can  be  concisely  written  as 

VCrt.pj.YpPyOa)  =  V(/?,PpYi  +~"*p2»a2)  ’• n  =  U.  (n> 

V(/?,PpYl,P2,a2)  =  VC^Pj.-Yj.pj-a,),  (12) 

VC^PpYi.PyOa)  =  V^PpYpTr-p^Tt+a,),  and  (13) 

V(R,pI,Y1.P2»a2)  =  V^pj.-YpTt-pj.Tt-^);  (14) 

further  symmetries  in  special  cases  can  be  written  as 

V(R,  PpYpO.o,)  =  V(R,  PpYj.0,0),  (15) 

V(i?,0,Y1,P2,a2)  =  V(/2,0,0,P2,a2-Yl),  (16) 

V(/?,7t,Y1,P2,a2)  =  V(i?,7t,0,p2,a2+Yi),  (17) 

VCR.O.YpO,^  =  V(i?,0,0,0,0),  and  (18) 

V(i?,7t,YpO,a2)  =  V(R,  7i,  0,0,0).  (19) 
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Including  these  symmetry  constraints  significantly  reduced  the  total  number  of  points  needed  for  our 
coarse  grid  covering  of  the  total  PES. 

A  total  of  187  points  on  the  PES  were  computed  for  various  combinations  of  /?,  P1?  yv  P2,  and 
Os-  The  first  group  of  points  chosen  was  based  on  the  four  geometries  presented  in  Table  2.  For  each 
of  these  geometries,  we  varied  only  the/?  coordinate  and  computed  enough  points  to  show  the  depth 
and  shape  of  the  local  minimum  well.  This  required  27  points  total  for  all  four  geometries 
investigated,  and  the  results  are  shown  in  Tables  6-8.  Next,  we  covered  a  coarse  grid  in  /?,  p,,  and 
Yj  while  fixing  the  remaining  coordinates  to  those  of  geometry  G1 .  For  this  collection  of  points,  we 
selected  values  for  R  from  the  set  (Gl,  3.5, 3.75, 4.0, 4.25, 4.5, 4.75  A).  For  px,  we  selected  from 
the  set  (0, 45, 90, 135, 180°)  and  for  yu  either  the  Gl  value  or  Gl  +60°.  A  subset  of  these  points, 
where  only  R  and  px  vary,  is  given  in  Tables  9  and  10  for  the  purpose  of  discussion,  and  the 
remainder  of  the  data  points  are  included  in  Tables  A12-A15  in  the  appendix.  A  clarification  is 
required  concerning  the  coordinates  for  the  points  in  Tables  A12-A15.  For  these  four  tables,  the 
algorithm  written  to  transform  from  internal  coordinates  to  Cartesian  coordinates  (used  as  input  to 
the  SAPT  procedure)  unintentionally  caused  a  reflection  through  the  XZ  plane  of  the  y1  value  in 
structure  Gl .  This  generated  Gl  -like  structures  that  differ  from  Gl  by  only  +0.294°  in  the  Yi  torsion 
angle.  Accordingly,  the  starting  yl  value  used  in  Tables  A12-A15  is  not  119.853°,  but  rather 
119.853°  +  0.294°,  or  taking  into  account  the  C3v  axis,  Yi  =  +0.147°.  The  coordinates  in 
Tables  A12-A15  include  this  difference  in  yv  Finally,  we  selected  a  coarse  grid  covering  the  entire 
PES  with  values  for/?  =  (3,4, 5, 6  A),  p!  =  (0, 45, 90, 135, 180°),  yi  =  (0, 60°),  p2  =  (0, 45°),  and 
ct2  =  (0, 45, 90°). 

Since  SAPT  interaction  energies  are  rigorously  free  of  BSSE,  there  is  no  a  priori  need  to 
compute  the  SAPT  components  in  a  dimer-centered  basis  set  (DCBS)  as  is  necessary  in  the  SM 
CP-coirected  approach.  Williams  et  al.  [38]  investigated  some  alternate  schemes  for  the  placement 
of  basis  functions  for  efficient  computation  of  SAPT  components,  called  monomer-centered  and 
monomer-centered  plus  basis  sets  (MCBS  and  MC+BS),  respectively.  In  a  “pure”  MCBS,  each 
monomer  uses  only  those  basis  functions  that  are  placed  on  its  own  nuclear  centers.  The  “plus”  in 
the  MC+BS  case  indicates  basis  functions  used  in  addition  to  the  monomers’  original  basis  set. 
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The  location  for  these  additional  functions  in  the  MC+BS  are  the  original  locations  of  the  basis 
functions  on  the  ghost  monomer.  Eventually,  adding  functions  to  the  MCBS  in  this  manner  would 
produce  the  full  DCBS.  The  goal,  however,  is  to  reduce  the  computational  effort  below  that  required 
for  the  full  DCBS  while  retaining  acceptable  accuracy.  A  MC+BS  that  best  balances  the  goals  of 
computational  tractability  and  accuracy  in  S  APT  calculations  is  one  in  which  only  the  valence  basis 
functions  are  retained  on  the  ghost  monomer.  Explicitly,  when  a  CH3CN,  SAPT  computation  is 
being  performed,  only  s-  and  p-type  functions  are  included  in  the  C02  ghost  basis  set.  Likewise, 
when  the  C02  computation  is  being  performed,  only  5  and  p  functions  are  placed  on  the  C  or  N 
atoms  of  CH3CN,  and  only  5-type  functions  are  placed  on  the  H  atoms.  This  arrangement  allows  a 
20-30%  reduction  in  the  size  of  the  basis  compared  to  the  equivalent  DCBS,  with  almost  no  sacrifice 
in  accuracy.  This  significantly  reduces  the  computational  effort  because  the  most  computationally 

expensive  SAPT  component,  E  ^ ,  scales  as  n30n*  where  and  nAv  are  the  numbers  of  occupied 

and  virtual  orbitals,  respectively.  In  this  manner,  the  full  DCBS  size  of  165  basis  functions  was 
reduced  to  a  MC+BS  size  of  135  and  117  basis  functions  for  monomer  CH3CN  and  C02, 
respectively. 

A  separate  computation  was  performed  for  the  component,  which  is  variational  in 

character.  The  MC+BS  was  augmented  with  a  large  midbond  set  of  2s2p2d\flg  placed  at  the 
midpoint  of  the  line  segment  defining  the  coordinate  R.  The  orbital  exponents  are  0.15  and  0.6  for 
5,  p,  and  d  functions  and  0.3  for  the/and  g  functions.  These  functions  were  selected  by  optimizing 
one  basis  function  of  each  type  to  one  digit  accuracy  with  the  goal  of  maximizing  the  absolute  value 

of  the  E  energy.  The  s,  p,  and  d  functions  were  then  each  split  into  two  basis  functions  using  the 

“even  scaling  rule”  [39]  producing  the  orbital  exponents  given  previously.  This  new  value  for  E  ^ 

was  then  used  in  place  of  the  one  which  was  computed  using  the  basis  set  described  in  the  previous 
paragraph.  The  addition  of  midbond  functions  for  this  component  will  be  indicated  with  a 
superscript  “( +mb)”  for  both  the  leading  term  in  the  second-order  dispersion  energy  and  for  the  total 

interaction  energy  as  £'^s°p(+mb)  and  respectively. 
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4.  Results  and  Discussion 


To  check  the  quality  of  the  atomic  orbital  (AO)  basis  set  used  on  the  monomers,  we  report  the 
calculated  electric  dipole  and  quadrupole  moment  of  COz  and  CH3CN,  respectively.  The  theoretical 
values  are  calculated  at  the  QCISD  level  with  the  aug-cc-pvdz  basis  set  used  on  the  monomers 
throughout  this  work.  The  electric  dipole  moment  of  CH3CN  is  predicted  to  be  3.95  debye,  while 
experimenting  gives  3.92  debye  [41  ] .  Values  for  the  C02  experimental  quadrupole  moment  [42-45] 
vary  between  - 1.34  and  1.5  x  10' 19  coulomb  A2,  while  the  calculated  value  lies  at  the  upper  end  of 
this  range  with  a  value  of  - 1.514  x  10' 19  coulomb  A2. 

Table  3  compares  the  S  APT  components  computed  using  the  MC  +BS  described  in  section  3  with 
the  equivalent  DCBS  at  the  two  geometries  G1  and  G4.  The  difference  in  the  first-order  components 

between  a  MC+BS  and  the  equivalent  DCBS  is  less  than  about  0.01  kcal/mol.  Aside  from  E^, 

which  will  be  replaced  as  previously  described,  the  second-order  components  differ  by  no  more  than 
0.04  kcal/mol.  Further,  the  second-order  differences  at  least  partially  offset  each  other  for  the 
geometries  shown,  though  this  cannot  be  relied  upon  over  the  entire  PES.  This  level  of  error  should 
be  substantially  smaller  than  other  sources  of  error  in  the  present  work.  Further,  the  reduction  in  the 
number  of  basis  functions  decreases  the  computational  cost  at  each  geometrical  configuration  by 

more  than  a  factor  of  three.  For  the  most  time-consuming  component,  the  triples  portion  of  E  ^ , 
this  cost  is  cut  by  3.5  times  by  using  a  MC+BS  rather  than  a  DCBS. 

Table  3  further  shows  that  the  largest  components  in  absolute  magnitude  for  G1  and  G4  are 
E  2® ,  e£°1,  and  EjJJj  •  Previous  experience  [38]  and  Table  3  indicate  that  the  first  two  components 

converge  very  quickly  with  basis  set  and  are  sufficiently  converged  with  the  current  MC +BS.  More 
attention  must  be  paid  to  E^  since  small  percentage  errors  in  this  component  can  translate  into 

relatively  large  absolute  errors  in  the  final  energies.  With  this  in  mind,  the  convergence  properties 
of  this  component  with  respect  to  basis  set  saturation  will  be  studied  in  more  detail  later. 
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Table  3.  Comparison  of  the  SAPT  Interaction  Energy  Components  for  the  CH3CN-C02 
Interaction  Using  a  MC+BS  and  the  Analogous  DCBS  at  the  Two  Local  Minimum 
Geometries  G1  and  G4  (The  Units  for  the  Energies  Are  kcal/mol.) 


]  Geometry 

G1 

G4 

Basis  Type 

MC+BS 

DCBS 

MC+BS 

DCBS 

Eao) 

-3.64 

-3.64 

-1.08 

-1.07 

e(1)  (3) 

C’elst,resp'^'/ 

0.53 

0.53 

0.11 

0.11 

Em 

3.88 

3.88 

1.09 

1.09 

0.25 

0.26 

0.20 

0.21 

f(20) 

"  disp 

-3.07 

-3.22 

-1.21 

-1.26 

-0.03 

-0.05 

-0.08 

-0.10 

^  disp 

-3.09 

-3.28 

-1.30 

-1.36 

e(20) 

ind,resp 

-2.01 

-2.05 

-0.26 

-0.26 

re- (22) 

•^ind 

0.05 

0.05 

-0.03 

-0.03 

p(20) 

^  exch -disp 

0.32 

0.36 

0.07 

0.08 

p(20) 

"  exch-ind,resp 

1.57 

1.57 

0.15 

0.15 

F  a 

^  int 

-2.38 

-2.51 

-1.06 

-1.12 

a  Computed  according  to  equation  (9). 


Table  4  compares  the  best  currently  coded  SAPT  approximation  to  each  SM  energy  (HF,  MP2, 
MP3,  and  MP4),  with  the  appropriate  sums  detailed  in  the  footnotes  to  this  table.  At  the  G1 

geometry  8  ,  the  difference  between  the  four-term  SAPT  approximation  to  Emt  and  E^  itself  is 

less  than  0.2  kcal/mol.  While  this  is  a  small  difference  in  absolute  magnitude,  it  represents  a  large 
percentage  of  E  ^ .  This  large  percentage  is  misleading,  though,  since  it  occurs  near  the  point  where  E 
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Table  4.  Comparison  of  SM  Results  With  S  APT  Results  at  Three  Local  Minimum  Geometries 
(The  SM  Interation  Energies  Are  Computed  Using  the  Boys-Bernardi  [25]  CP 
Scheme;  the  SAPT  Components  Use  a  MC+BS.  The  DCBS  Values  Available  for 
Geometries  G1  and  G4  Are  Displayed  in  the  Footnotes  if  They  Differ  by  More  Than 
0.02  kcal/mol  From  the  MC  +BS  Value.) 


Geometry 

Basis  Type 

G1 

G2 

G4 

£HF 
"  rnt 

-0.39 

1.30 

-0.15 

77  a 

"  SAPT 

-0.21 

1.36 

-0.10 

gBff 

-0.18 

-0.06 

-0.06 

Z?MP2 

^mt 

-1.75 

-1.49 

-0.75 

27  MP2  b 
"  SAPT 

-1.73 c 

-1.52 

-0.80 d 

gMP2 

-0.03 

0.05 

0.05 

Z7  MP3 

^  int 

0.47 

0.23 

0.09 

-MP3  e 
"  SAPT 

0.56 

0.30 

0.17 

gMP3 

-0.08 

-0.07 

-0.08 

fMP4 

int 

-0.50 

-0.30 

-0.16 

PMP4  f 
"  SAPT 

-0.81  * 

-0.30 

-0.29 

gMP4 

0.30 

0.00 

0.13 

r,  (HF  +  MP2  +  MP3  +  MP4) 

^  int 

-2.17 

-0.26 

-0.97 

E  h 

^  SAPT 

-2.37 

-0.22 

-1.07 

g  (HF  +  MP2  +  MP3  +  MP4) 

0.20  ; 

-0.04 

0.10 

a  E.HF  _  p(10)  .  P(10)  .  p.(20)  p(20) 

^  SAPT  ^  eist  ^  exch  ^  ind,resp  ^  exch-ind,resp  * 

b  -MP2  _  -(20)  p(12)  tp (22)  (1)  n v  17(20) 

^  SAPT  ^  disp  ^elsLresp  ^  ind.resp  ^exch^'  ^exch-disp 

c  DCBS  value  is  - 1.84  kcal/mol. 
d  DCBS  value  is  -0.84  kcal/mol. 

erMP3  _  p(  13)  77(21) 

^  SAPT  ~  "  elst.resp  ^disp* 

f  17MP4  _  17  (22) 

SAPT  ~  ^  disp’ 

8  DCBS  value  is  -0.83  kcal/mol. 
h  Computed  according  to  equation  (9). 


+ 


tp  (22) 

^  exch-ind* 
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crosses  zero.  The  absolute  values  of  and  E^h  given  in  Table  3  are  nearly  4  kcal/mol  at  Gl, 

but  are  opposite  in  sign  and  cancel  to  within  a  fraction  of  a  kcal/mol.  The  comparison  to  SM-MP2 
is  quite  good  at  all  three  geometries,  probably  indicating  reasonably  good  agreement  over  the  entire 
PES.  The  SM-MP3  energies  show  a  small  difference  in  absolute  values  of  less  than  0.1  kcal/mol. 
The  SM-MP4  energy  comparison  shows  the  largest  deviation  of  0.3  kcal/mol  from  the  SAPT 
approximation  at  geometry  Gl.  This  deviation  is  probably  due  to  accumulation  of  errors  resulting 

from  the  neglect  of  E<|Jresp  and  E^_disp. 


Table  5  investigates  the  effect  on  Ed-S°p  of  adding  basis  functions  with  different  angular 

symmetries  placed  at  the  midbond  position  as  described  in  section  3.  The  value  for  this  component 
using  the  DCBS  from  Table  3  is  repeated  in  Table  5  for  comparison.  The  addition  of  the  full 
midbond  2s2p2d\flg  set  lowers  this  energy  component  by  0.53  (17%)  and  0.18  (15%)  kcal/mol  for 
the  geometries  Gl  and  G4,  respectively,  from  its  MC+BS  value.  This  results  in  a  22%  and  16% 
lowering  of  the  final  value  of  E-mx  for  these  geometries.  By  adding  only  s  and  p  functions  at  the 

midbond,  the  value  of  E  ^s°p  is  already  lower  than  the  corresponding  DCBS  value  while  using  22 

fewer  basis  functions.  This  indicates  that  midbond  functions  are  helping  to  converge  this  component 
faster  than  nuclear-centered  basis  functions  of  higher  angular  symmetry  placed  at  the  other  nuclear 

centers.  The  increase  in  the  computational  cost  of  treating  only  E^s°p  with  the  full  set  of  midbond 

basis  functions  is  rather  minor,  yet  provides  a  significant  increase  in  accuracy  for  this  component. 
Using  this  set  of  midbond  functions  in  the  computation  of  all  SAPT  components  would  not  have 
improved  the  accuracy  of  the  results  in  proportion  to  the  amount  of  additional  computer  time  needed. 

Figure  3  displays  a  cut  in  R  for  each  of  the  four  geometries,  Gl,  G2,  G3,  and  G4.  Tables  6, 7, 
and  8  provide  the  individual  SAPT  components  and  total  interaction  energy  for  each  of  the  points 
shown  in  Figure  3.  Geometry  Gl ,  with  C02  roughly  parallel  to  the  CCN  axis,  has  a  total  interaction 
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Table  5.  Contribution  of  Basis  Functions  of  Different  Angular  Symmetry  Placed  at  the 
Midbond  Position  to  the  E  Energy  in  kcal/mol  (The  Table  Indicates  the  Successive 

Inclusion  of  Each  of  the  Symmetries,  e.g.,  Row  g  Includes  the  Full  2s2p2dlflg  Set. 
The  Numerical  Values  of  the  Exponents  Used  for  Each  Symmetry  Are  Found  in 
Section  3.  The  Number  of  Basis  Functions  Refers  to  the  CH3CN  Monomer  Including 
Midbond  Functions  and  Ghost  Functions  Placed  on  the  C02  Nuclear  Centers.  The 
Values  Computed  Without  Midbond  Functions  Are  Taken  From  Table  3.) 


Basis  Type 

Midbond 

Basis  Size 

Geometry  jj 

Number 

Symmetry 

G1 

G4 

DCBS 

— 

— 

— 

165 

-3.22 

-1.26 

MC+BS 

— 

— 

— 

135 

-3.07 

-1.21 

MC+BS 

2 

— 

s 

137 

-3.17 

-1.23 

MC+BS 

2 

— 

p 

143 

-3.35 

-1.31 

MC+BS 

2 

— 

d 

153 

-3.47 

-1.36 

MC+BS 

1 

— 

f 

160 

-3.55 

-1.38 

MC+BS 

1 

— 

§ 

169 

-3.60 

-1.39 

energy  of  Etot(+I“b)  =  -2.90  kcal/mol.  This  is  slightly  more  stable  than  G3,  where  Eint(+mb)  = 

-2.82  kcal/  mol,  and  the  C02  is  perpendicular  to  the  CCN  axis.  The  interaction  energy  of 
-1.25  kcal/mol  at  geometry  G4  is  significantly  smaller  in  magnitude  than  for  the  former  two 
geometries.  Finally,  geometry  G2  is  only  weakly  bound  at  -0.43  kcal/mol. 

Figure  4a  shows  the  cut  through  geometry  G1  broken  into  the  components  indicated  in 
equation  (9).  The  minimum  energy  for  E^f  is  near  R  =  3.75  A,  while  the  full  interaction  energy, 

including  correlation  corrections,  predicts  the  minimum  to  occur  near  R  =  3.33  A  or  0.4  A  shorter 
than  the  HF  value.  This  comparison  clearly  indicates  the  need  to  include  the  intermolecular  electron 
correlation  effects  for  this  system.  At  the  minimum  for  the  correlated  interaction,  the  HF 
contribution  is  only  -0.39  kcal/mol  and  is  determined  by  two  pairs  of  large  values  of  opposite  signs: 

the  first-order  interactions  E^j,  +  £ dS?  =  3.88  -  3.62  =  +0.26  kcal/mol  and  the  second-order 
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Figure  3.  Cuts  in  R  Through  the  PES  for  the  Four  Local  Minimum  Geometries  Detailed  in  Table  2.  The  Final  Energies  Were 
Computed  Using  Equation  (9)  and  Include  E  2sp(+mb)  •  The  Solid  Line  Is  a  Spline  Fit  to  the  Single-Point  SAPT  Energies 

for  Each  Cut  and  Is  Only  to  Guide  the  Eye.  Energies  Are  in  kcal/mol,  and  Distances  Refer  to  the  Center-of-Mass 
Separation  Between  the  Two  Monomers  in  A. 


Table  6.  Potential  Cuts  in  R  Through  the  Two  Local  Minimum  Geometries  G1  and  G2  (All  Energies  Are  in  kcal/mol  and  Distances 
in  A.) 
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Table  7.  Potential  Cuts  in  R  Through  the  Local  Minimum  Geometry  G3  (All  Energies  Are  in 
kcal/mol  and  Distances  in  A.) 


1  R  |  3.75  |  4.00 

1  4.28  1  4.50  1  4.75 

5.00  |  5.50  | 

E™ 

int 

5.03 

0.54 

-1.13 

-1.44 

-1.39 

-1.20 

-0.82 

77  (10) 

"  elst 

-11.29 

-6.16 

-3.52 

-2.43 

-1.72 

-1.28 

-0.79 

r(10) 

^  exch 

18.84 

7.87 

2.96 

1.32 

0.53 

0.21 

0.03 

77  (20) 

"  ind.resp 

-9.30 

-3.71 

-1.42 

-0.69 

-0.33 

-0.18 

-0.07 

EW) 

exch -ind,  resp 

7.48 

2.84 

0.96 

0.40 

0.14 

0.05 

0.01 

gHF 

-0.70 

-0.29 

-0.10 

-0.04 

-0.02 

-0.01 

0.00 

£0)  (2) 
^elst.resp'^' 

-0.30 

0.01 

0.15 

0.18 

0.17 

0.16 

0.11 

02) 

2.13 

1.14 

0.53 

0.28 

0.13 

0.06 

0.01 

tF  (22) 

"  ind 

-1.17 

-0.51 

-0.19 

-0.08 

-0.03 

-0.01 

0.00 

lr  (22) 

^  exch -ind 

0.94 

0.39 

0.13 

0.05 

0.01 

0.00 

0.00 

17(20) 
c  disp 

-5.91 

-3.55 

-2.06 

-1.36 

-0.87 

-0.58 

-0.28 

17  (21) 

"  disp 

1.51 

0.87 

0.50 

0.33 

0.22 

0.15 

0.07 

17(22) 

"  disp 

-1.72 

-1.03 

-0.61 

-0.40 

-0.26 

-0.17 

-0.08 

eSp(2) 

-0.21 

-0.16 

-0.11 

-0.07 

-0.04 

-0.02 

-0.01 

-6.12 

-3.71 

-2.17 

-1.43 

-0.92 

-0.60 

-0.28 

e(20) 
exch -disp 

1.15 

0.54 

0.23 

0.11 

0.04 

0.02 

0.00 

-3.37 

-2.14 

-1.33 

-0.90 

-0.58 

-0.38 

-0.15 

^int 

1.66 

-1.60 

-2.45 

-2.34 

-1.97 

-1.58 

-0.97 

77  (20)  (+ mb) 

^  disp 

-7.13 

-4.24 

-2.43 

-1.58 

-1.00 

-0.66 

-0.30 

17  (-mb) 

int 

0.44 

-2.29 

-2.82 

-2.57 

-2.10 

-1.66 

-1.00 
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Table  8.  Potential  Cuts  in  R  Through  the  Local  Minimum  Geometry  G4  (All  Energies  Are  in 
kcal/mol  and  Distances  in  A.) 


R 

4.75 

5.00 

5.25 

5.57 

5.75 

6.00 

6.50 

rHF 

E  int 

14.29 

5.11 

1.43 

-0.15 

-0.43 

-0.54 

-0.48 

Em 

c  elst 

-8.55 

-3.89 

-2.00 

-1.08 

-0.85 

-0.66 

-0.46 

■H 

24.74 

9.84 

3.82 

1.09 

0.54 

0.20 

0.03 

p  (20) 

^  ind,resp 

-6.43 

-2.22 

-0.81 

-0.26 

-0.15 

-0.08 

-0.04 

Fm 

£''excb-md,resp 

5.67 

1.92 

0.64 

0.15 

0.07 

0.02 

0.00 

5^ 

-1.15 

-0.53 

-0.21 

-0.05 

-0.02 

-0.01 

0.00 

e(1)  (3) 

fcelst,resp 

-0.48 

-0.12 

0.04 

0.11 

0.12 

0.12 

0.09 

«2»<2) 

2.16 

1.15 

0.56 

0.20 

0.11 

0.05 

0.01 

r*r(  22) 

^ind 

-0.99 

-0.38 

-0.14 

-0.03 

-0.01 

0.00 

0.00 

/F(  22) 
c  exch-ind 

0.88 

0.33 

0.11 

0.02 

0.01 

0.00 

0.00 

,7  (20) 

c  disp 

-6.38 

-3.74 

-2.25 

-1.21 

-0.88 

-0.58 

-0.27 

77(21) 

disp 

1.19 

0.65 

0.38 

0.21 

0.16 

0.11 

0.05 

e(22) 

^  disp 

-1.30 

-0.80 

-0.50 

-0.29 

-0.22 

-0.15 

-0.08 

eSP(2) 

-0.11 

-0.14 

-0.12 

-0.08 

-0.06 

-0.04 

-0.02 

*Sp© 

-6.49 

-3.88 

-2.37 

-1.30 

-0.95 

-0.62 

-0.29 

Em 

^  exch-disp 

1.12 

0.51 

0.22 

0.07 

0.04 

0.62 

0.00 

-3.80 

-2.39 

-1.57 

-0.92 

-0.68 

-0.44 

-0.19 

*in, 

10.49 

2.73 

-0.14 

-1.07 

-1.11 

-0.98 

-0.66 

p  (20)(+mb) 

"  disp 

-7.94 

-4.55 

-2.66 

-1.39 

-0.99 

-0.64 

-0.29 

17  (+mb) 

"  int 

8.93 

1.92 

-0.55 

-1.25 

-1.22 

-1.04 

-0.68 
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£exch-ind,resp  +  ^^d!resp =  1-57  -  2.01  =  -0.44kcal/mol.  The  resulting  two  values  of  opposite  signs, 

plus  the  residual  e  =  -0.18  [see  equation  (8)],  produce  the  comparatively  small  -0.39  kcal/mol 
contribution  by  the  Em .  Therefore,  the  correlated  interaction  terms  account  for  most  of  the 

-2.90  kcal/mol  stabilization  energy  £'int(+mb)  near  the  minimum.  It  is  worth  pointing  out  that  even 

though  the  net  HF  induction  energy  is  only  -0.44  kcal/mol,  it  still  represents  15%  of  the  total 
£int<+mb)  =  -2.90  kcal/mol. 

Closer  examination  of  Figure  4a  shows  that  almost  all  of  the  stabilizing  energy  comes  from  the 
dispersion  energy.  The  intramonomer  electron  correlation  contributions,  and  E  ^sp ,  essentially 

cancel  one  another,  making  a  very  small  net  attractive  contribution  [the  sum  is  labeled  as 

Table  6  shows  that  at  the  G1  conformation,  €(^sp(2)  contributes  only  1%  of  the  total  dispersion 

energy.  In  general,  €®sp(2)  accounts  for  no  more  than  about  10%  of  the  leading-order  component, 

E a°p(+mb) ,  and  for  about  two-thirds  of  all  of  the  points  computed,  this  contribution  is  less  than  5%. 

The  primary  source  of  the  stabilization  of  the  G1  complex  can  be  pinpointed  to  the  leading-order 
dispersion  term  E  mb) . 

Moving  away  from  the  G1  minimum  to  larger  intermolecular  distances,  the  stabilization  energy 
at  R  =  3.75  A  (see  Figure  4a)  is  composed  of  nearly  equal  contributions  from  the  dispersion  term  and 

the  HF  energy.  The  delicate  balance  between  the  first-  ^  and  E^l)  and  second-  (Zij^resp  and 
£’exch-ind.resp) order  interactions  contributing  to  is  shown  in  Figure  4b.  Table  6  shows  that  the 
repulsive  first-order  exchange  energy  is  a  major  contributor  to  over  most  of  the  cut  through  G1 
and  counters  the  large  stabilizing  interaction  arising  from  the  electrostatic  term  E  ^ .  This  gives  a 
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E  [kcal/mol] 


Figure  4a.  SAPT  Components  Used  in  Equation  (9)  for  a  Range  of  R  Values  at  and  Around 
the  Local  Minimum  Geometry  Gl.  The  Estimated  Component  {E ^ch-ind  Is  Not 

Included  in  the  Figure.  The  Numerical  Results  for  This  Cut  Are  Displayed  in 
Table  6.  The  Energies  and  Distances  Are  Given  in  Units  of  kcal/mol  and  A, 
Respectively. 


25 


E  [kcal/mol] 


Figure  4b. 


E  and  the  SAPT  Components  Included  in  Equation  (8)  at  and  Around 
Geometry  G1  as  a  Function  of  R. 


26 


net  HF  first-order  contribution  of  =  -0.93  kcal/mole,  or  39%  of  the  total  interaction 

energy  at  R  =  3.75  A.  While  this  is  also  a  general  trend,  the  G2  geometry  is  a  notable  exception. 
Table  6  shows  the  G2  HF  electrostatic  term  E  ^  to  be  repulsive,  with  a  value  of  +0.33  kcal/mol 

near  the  minimum  at  R  =  4.65  A.  Returning  to  the  G1  geometry,  the  induction  energy  (F^.resp) 

accounts  for  22%  of  the  stabilization  interactions  ati?  =3.75  A,  which  gives  a  net  HF  second-order 
contribution  of  F^resp  +  ^ixch-ind.resp  =  -0.19,  or  8%  of  the  total  interaction  energy. 

A  summary  of  the  individual  components  split  according  to  equation  (9)  for  all  four  geometries, 
G1 ,  G2,  G3,  and  G4,  are  shown  in  histogram  format  in  Figure  5a.  A  quick  inspection  shows  that  the 

is  disp  (2) ( + mb)  energy  is  the  largest  component  in  absolute  value  for  each  of  these  configurations. 
Also,  it  is  clear  that  plays  an  important  role  in  determining  the  total  interaction  energy  at  two 
of  the  four  minimum  geometries,  specifically  G2  and  G3.  For  all  configurations  except  G2,  the  E^ 
energy  stabilizes  the  complex. 

HF 

A  histogram  for  the  components  of  for  these  four  geometries  is  given  in  Figure  5b.  The 
repulsive  contribution  of  for  G2  differs  from  the  other  three  minima  by  having  a  positive 
electrostatic  interaction  energy  (i.e.,  E^  =  +0.33  kcal/mol).  In  the  other  three  minimum 
conformations,  is  negative  and  nearly  cancels  the  positive  F^ch  term.  Hence,  as  seen  in 

Table  6,  the  combination  of  the  positive  electrostatic  term  and  a  large  (positive)  E^h  account  for 

the  very  weak  bond  at  this  geometry.  Finally,  for  configurations  G1  and  G4,  Figure  5b  illustrates 
the  cancellation  within  the  first-order  and  second-order  HF  terms,  leaving  02)(+mb)  as  Ae 

dominant  contribution  to  the  total  interaction  energy. 
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Figure  5a.  The  Interaction  Energy  and  Individual  Components  Indicated  in  Equation  (9)  Are  Shown  in  Histogram  Format  to 
Illustrate  the  Relative  Importance  of  the  Various  Energy  Contributions  at  the  Gl,  G2,  G3,  and  G4  Minima. 
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Figure  5b.  A  Histogram  Is  Shown  Coni 
G4  Minima. 


Tables  9  and  10  provide  the  interaction  energies  on  the  PES  around  geometry  G1  as  a  function 
of  Pi  for  three  values  of  R.  Figure  6  displays  the  p3  dependence  of  the  total  interaction  energy  at 
R  =  4.0  and  4.75  A.  Figure  6  clearly  shows  the  asymmetry  in  Pj  for  the  R  =  4.0  A  cut  due  to  the 
absence  of  a  oh  symmetry  plane  in  CH3CN.  At  R  =  4.75  A,  the  PES  flattens  out  and  is  attractive 
everywhere  as  a  function  of  the  single  variable  Pt,  albeit  weakly  bound.  These  plots  predict  a 
favorable  angle  P!  of  approach  in  the  range  80°  ^  Pi  <  150°. 

Analyses  of  the  terms  defining  the  total  Fint(+mb)  are  shown  in  Figures  7  and  8  for  R  =  4.0  and 
4.75  A,  respectively.  Figure  7a  plots  the  components  of  [according  to  equation  (9)]  as  a 

function  of  Pj  at  R  =  4.0  A.  Again,  Fint("'mb)  is  determined  primarily  by  two  terms,  F*f  and 
F^p(2)(+mb).  The  remaining  terms  from  Ec™  are  small  in  absolute  magnitude  and  essentially 
cancel  one  another.  Table  9  indicates  that  the  leading-order  dispersion  term,  F^sp(+mb),  accounts 
for  most  of  the  net  contribution  to  the  total  dispersion  interaction.  The  E  ^p  and  E  ^  energies  are 

individually  28%  or  less  of  the  F^°p  energy,  but  are  nearly  equal  in  magnitude  and  opposite  in  sign, 

effectively  canceling  one  another.  This  cancellation  of  dispersion  components  was  noted  earlier  in 
the  analysis  of  Fint(+mb)  as  a  function  of  R. 

The  separation  of  F®  into  the  components  indicated  in  equation  (8)  is  shown  in  Figure  7b.  The 

second-order  terms,  F^.^  resp  and  F^resp,  practically  cancel  over  the  entire  range  of  Pt.  Thus, 

the  general  dependence  of  F^f  on  Pj  can  be  assigned  to  the  two  first-order  terms  F^  and  E^ . 

As  Pj  goes  to  180°,  the  atoms  of  the  C02  molecule  come  into  closer  proximity  with  atoms  on 
CH3CN.  At  Pj  =  180°,  an  oxygen  resides  only  1 .95  A  away  from  the  nitrogen.  In  this  arrangement, 
one  would  expect  a  repulsive  intermolecular  penetration  of  the  electron  clouds  by  the  two  monomers 
(i.e.,  a  “steric”  interaction,  and  indeed  the  exchange  [repulsion]  interaction  dominates  near  180°). 


30 


Table  9.  Potential  Cuts  in  Pj  at  Geometry  G1  for  Two  Different  R  Values  (All  Energies  Are 
in  kcal/mol  and  Distances  in  A.  Only  Two  Angles  Are  Shown  for  thei?  Value  Taken 
From  Geometry  Gl,  Since  Each  of  the  Others  Is  Too  High  on  the  Exponential  Wall.) 


R 

Gl 

4.0 

Eli 

90° 

135° 

45° 

90° 

108° 

135° 

180° 

rHF 

E-m 

1.58 

4.16 

11.04 

-0.8 

-1.14 

-1.01 

9.59 

E(10) 

celst 

-3.5 

-8.45 

-6.13 

-1.14 

-1.40 

-2.41 

-8.92 

rp 

^  exch 

5.82 

14.31 

18.46 

0.5 

0.42 

1.74 

20.09 

Fm 

c  ind.resp 

-1.91 

-7.76 

-5.64 

-0.23 

-0.27 

-0.9 

-8.22 

e(20) 

exch -ind,  resp 

1.50 

6.58 

5.16 

0.1 

0.14 

0.64 

7.24 

gHF 

-0.32 

-0.53 

-0.81 

-0.03 

-0.03 

-0.08 

-0.59 

e(1)  (3) 

c  elst,resp 

0.24 

0.43 

-0.32 

0.21 

0.29 

0.29 

-0.56 

0.58 

0.68 

1.13 

0.08 

0.03 

0.13 

1.52 

ip  (22) 

^  ind 

-0.14 

-0.06 

-0.64 

0.00 

0.02 

0.01 

-1.26 

•e(22) 

exch  -ind 

0.11 

0.05 

0.58 

0.00 

-0.01 

-0.01 

1.11 

e(20) 

^  disp 

-3.57 

-5.55 

-5.13 

-1.09 

-0.98 

-1.64 

-5.57 

e(21) 

^  disp 

0.80 

1.48 

0.78 

0.25 

0.25 

0.41 

1.46 

e{22) 

^  disp 

-0.88 

-1.43 

-0.96 

-0.29 

-0.27 

-0.42 

-1.55 

4P(2) 

-0.08 

-0.05 

-0.17 

-0.03 

-0.01 

-0.01 

-0.1 

O2) 

-3.65 

-5.5 

-5.3 

-1.12 

-0.99 

-1.65 

-5.66 

Ee°) 

exch  -disp 

0.40 

1.01 

0.95 

0.04 

0.04 

0.15 

1.22 

** 

-2.46 

-3.38 

-3.61 

-0.79 

-0.62 

-1.08 

-3.62 

^int 

-0.88 

0.78 

7.43 

-1.59 

-1.76 

-2.09 

5.97 

rp  (^O)  (  +  mb) 

^  disp 

-4.11 

-6.74 

-6.04 

-1.21 

-1.12 

-1.95 

-6.75 

t?  (+mb) 

^  int 

-1.41 

-0.41 

6.52 

-1.71 

-1.89 

4.79 
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Table  10.  Potential  Cuts  in  pt  at  Geometry  G1  for  R  =  4.75  A  (All  Energies  Are  in  kcal/mol.) 
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Figure  7a.  The  Interaction  Energy  E  [nt"'b)  and  Its  Components  According  to  Equation  (9)  Shown  as  a  Function  of  fJ 
The  Remaining  Coordinates  Are  Those  of  the  G1  Geometry.  (Data  Taken  From  Table  9.) 
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Figure  7b.  The  HF  Interaction  Energy  E^x  and  Its  Components  According  to  Equation  (8)  Shown  as  a  Function  of  px  for 
R  =  4.0  A.  The  Remaining  Coordinates  Are  Those  of  the  G1  Geometry.  (Data  Taken  From  Table  9.) 
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The  Interaction  Energy  E\ni  and  Its  Components  According  to  Equation  (9)  Shown  as  a  Function 
R  =4.75  A.  The  Remaining  Coordinates  Are  Those  of  the  G1  Geometry.  (Data  Taken  From  Table  10.) 
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Figure  8b.  The  HF  Interaction  Energy  E  and  Its  Components  According  to  Equation  (8)  Shown  as  a  Function  of  P,  for 
R  =  4.75  A.  The  Remaining  Coordinates  Are  Those  of  the  G1  Geometry.  (Data  Taken  From  Table  10.) 


This  strong  dependence  of  on  (indirectly  on  R )  can  be  seen  in  Tables  9  and  10  and  in 

Figure  7b  as  changes  in  Pj  bring  the  two  molecules  closer  together.  The  *erm  charges  from 
a  maximum  of  +20  kcal/mol  at  px  =  180°  to  0.4  kcal/mol  at  108°  (see  Table  9). 

As  the  two  molecules  separate  to  a  distance  of  R  =  4.75  A,  Figure  8a  shows  a  similar  relationship 
between  the  total  interaction  energy  and  its  components  as  seen  at  the  shorter  R  =  4.0  A  distance. 

As  Pj  varies  over  45°  ^  Pj  ^  180°,  £’int(+mb)  is  once  again  seen  to  be  composed  primarily  of  the 
and  E®  (2)<+mb>  contributions,  and  the  dispersion  interaction  is  comprised  almost  entirely  of  the 
leading-order  term  £^°p(+mb)  (see  Table  10).  E^f  is  again  a  balance  between  the  first-order  terms  E  ^ 

and  E  ^ .  From  the  previous  discussion,  it  is  clear  that  the  qualitative,  if  not  quantitative,  changes 
in  E  int(+mb)  as  a  function  of  R  and  px  can  be  traced  primarily  to  three  interaction  terms:  E  ^ ,  E  ^ , 
and  2?g(+mb>. 

5.  Conclusions 

The  CH3CN-C02  PES  was  investigated  using  SAPT.  Approximately  200  geometrical 
configurations  were  computed  for  both  selected  individual  cuts  of  the  PES  around  the  4  minimum 
energy  geometries  and  a  coarse  grid  spanning  the  5  intermolecular  coordinates.  A  separate 
computation  with  a  larger  basis  set,  including  midbond  functions,  was  also  performed  at  each  point 
for  the  leading-order  dispersion  energy.  Four  representative  local  minimum  geometries  were 
investigated  to  determine  the  relative  strengths  of  different  physical  contributions  to  the  interaction 

energy.  The  leading-order  dispersion  energy,  E*^,  contributed  a  large  percentage  of  the  binding 

energy  near  each  of  the  local  minima  investigated.  Surprisingly,  the  intramonomer  electron 
correlation  corrections  to  the  leading-order  dispersion  component  had  very  little  impact  on  the  final 
energies  due  to  cancellation  between  them — even  though,  individually,  they  typically  had  large 


absolute  magnitudes.  The  SM  HF  energy,  which  includes  contributions  from  the  leading-order 

SAPT  components  of  polarization,  exchange,  and  induction,  also  had  a  large  though  varying  impact 

HF 

on  the  final  interaction  energy  around  the  local  minima  investigated.  The  main  contributors  to  E  ^ 

are  the  first-order  exchange  and  electrostatic  terms,  where  the  importance  of  the  electrostatic  term 
is  due  to  the  strongly  polar  nature  of  the  CH3CN  monomer.  The  most  strongly  bound  geometrical 
configuration  investigated  had  an  interaction  energy  of  -2.90  kcal/mol. 
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Appendix: 

Supporting  Information  for  Investigation  of  the 
CH3CN-C02  Potential  Energy  Surface  (PES)  Using 
Symmetry- Adapted  Perturbation  Theory  (SAPT) 
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Table  A-l.  The  Supplementary  Tables  Complement  Tables  6-10  and  Provide  a  Complete  Description  of  All  Single-Point 

Symmetry-Adapter  Perturbation  Theory  (SAPT)  Computations  Used  in  the  Present  Work  (The  Units  for  Energies, 
Distances,  and  Angles  Are  kcal/mol,  A,  and  Degrees,  Respectively.) 
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Table  A-l.  The  Supplementary  Tables  Complement  Tables  6-10  and  Provide  a  Complete  Description  of  All  Single-Point 

Symmetry-Adapter  Perturbation  Theory  (SAPT)  Computations  Used  in  the  Present  Work  (The  Units  for  Energies, 
Distances,  and  Angles  Are  kcal/mol,  A,  and  Degrees,  Respectively.)  (continued) 
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Table  A-2.  Supplementary  Tables  Continued  (For  Description  and  Units,  See  the  Caption  of  Table  A-l.) 
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Supplementary  Tables  Continued  (For  Description  and  Units,  See  the  Caption  of  Table  A-l.)  (continued) 
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Table  A-8.  Supplementary  Tables  Continued  (For  Description  and  Units,  See  the  Caption  of  Table  A-l.) 
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.  Supplementary  Tables  Continued  (For  Description  and  Units,  See  the  Caption  of  Table  A-l.)  (continued) 
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Table  A-12.  Supplementary  Tables  Continued  (For  Description  and  Units,  See  the  Caption  of  Table  A-l.) 
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Full  coordinate  value  is  60.147002. 
Full  coordinate  value  is  0.147002. 


Table  A-12.  Supplementary  Tables  Continued  (For  Description  and  Units,  See  the  Caption  of  Table  A-l.)  (continued) 
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Table  A-13.  Supplementary  Tables  Continued  (For  Description  and  Units,  See  the  Caption  of  Table  A-l;  See  Footnotes  in 
Table  A-12  for  Coordinate  Information.) 
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